Mineria de dades: Earthquakes Part 2
- 1 Introducció
- 2 Exercici 1
- 3 Exercici 2
- 4 Exercici 3
- 5 Exercici 4
- 6 Exercici 5
- 6.1 Creació del model i visualització de l’arbre.
- 6.2 Extracció de les regles, comentaris i interpretacions.
- 6.3 Qualitat del model
- 6.4 Generació la matriu de confusió per a mesurar la capacitat predictiva de l’algoritme, tenint en compte les diferents mètriques associades a aquesta matriu (precisió, sensibilitat, especificitat).
- 6.5 Conclusions.
- 7 Exercici 6
- 8 Exercici 7
1 Introducció
Com a continuació de l’estudi iniciat a la Part 1, procedim a aplicar models analítics, tant no supervisats com supervisats, sobre el joc de dades seleccionat i preparat. En aquesta Part 2 partim de les dades prèviament preparades a la Part 1.
Models no supervisats
Aplicar un model no supervisat basat en el concepte de distància, sobre el joc de dades.
Aplicar de nou el model anterior, però usant una mètrica de distància diferent i comparar-ne els resultats amb els mètodes anteriors.
Utilitzar els algorismes DBSCAN i OPTICS, provant amb diferents valors del paràmetre
epsiminPts, i es comparen els resultats amb els mètodes anteriors.
Models supervisats
Seleccionar una mostra d’entrenament i una de test utilitzant les proporcions que es considerin més adequades en funció de la disponibilitat de dades. Justificar aquesta selecció.
Aplicar un model de generació de regles a partir d’ arbres de decisió ajustant les diferents opcions de creació. Obtenir l’arbre sense i amb opcions de poda. Obtenir la matriu de confusió. Finalment, comparar-ne els resultats.
Aplicar un model supervisat diferent del del punt 5., s’ha de triar entre els que s’han vist al material docent de l’assignatura. Comparar el resultat amb el model generat anteriorment.
Identificar eventuals limitacions del dataset seleccionat i analitzar els riscos per al cas d’ús del model per a classificar una nova dada.
1.1 Conjunt de dades
Abans de començar amb els exercicis, recordem quines són les nostres dades i en quin estat les tenim.
El nostre joc de dades procedeix de Kaggle, en particular:
Global Significant Earthquake Database from 2150BC, és un llistat d’uns 6000 terratrèmols que van del 2150 BC fins a l’actualitat. Ja vam veure que aquest conjunt de dades és altament incomplet. Un cop hem netejat els valors NA’s tenim:
library(tidyverse)
dades <- read_csv('Worldwide-Earthquake-database.csv')# Fem la neteja de la pràctica 1
dades <- dades |>
select(1:31)
# Considerem terratrèmols del 1900 endavant
dades_1900 <- dades |>
filter(YEAR >= 1900)
# Eliminem NA's de la profunditat i magnitud
# Per Pra2 no eliminem valors NA's de hora i morts, com vam fer a la Pra1
dades_netes <- dades_1900 |>
filter(!is.na(FOCAL_DEPTH) & !is.na(EQ_PRIMARY) & !is.na(DAMAGE_DESCRIPTION) & !is.na(DEATHS))
# Convertim Yes = 1, No = 0
dades_netes <- dades_netes |>
mutate(FLAG_TSUNAMI = ifelse(FLAG_TSUNAMI %in% c("No"), 0, 1))Per tant, de moment ens quedem amb:
# variables amb les que ens quedem (2 = FLAG_TSUNAMI, 18 = COUNTRY, no la posem, so far)
# cols <- c(2, 3, 9, 10, 18, 21, 22, 24, 25, 31)
# dades_countries <- as.data.frame(dades_netes[cols])
cols <- c(2, 3, 9, 10, 21, 22, 24, 25, 31)
dades <- as.data.frame(dades_netes[cols])library(skimr)
skim_without_charts(dades)| Name | dades |
| Number of rows | 1137 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| numeric | 9 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 |
|---|---|---|---|---|---|---|---|---|---|
| FLAG_TSUNAMI | 0 | 1 | 0.21 | 0.41 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 |
| YEAR | 0 | 1 | 1984.57 | 27.55 | 1900.00 | 1969.00 | 1990.00 | 2006.00 | 2020.00 |
| FOCAL_DEPTH | 0 | 1 | 31.63 | 40.53 | 0.00 | 10.00 | 22.00 | 33.00 | 651.00 |
| EQ_PRIMARY | 0 | 1 | 6.38 | 0.97 | 2.10 | 5.70 | 6.40 | 7.10 | 9.50 |
| LATITUDE | 0 | 1 | 19.71 | 21.46 | -54.00 | 2.71 | 27.50 | 37.19 | 61.02 |
| LONGITUDE | 0 | 1 | 38.36 | 79.55 | -178.25 | -0.66 | 52.19 | 103.32 | 178.29 |
| DEATHS | 0 | 1 | 1872.21 | 15359.15 | 1.00 | 2.00 | 8.00 | 60.00 | 316000.00 |
| DEATHS_DESCRIPTION | 0 | 1 | 1.56 | 1.01 | 1.00 | 1.00 | 1.00 | 2.00 | 4.00 |
| DAMAGE_DESCRIPTION | 0 | 1 | 2.51 | 1.04 | 1.00 | 2.00 | 2.00 | 3.00 | 4.00 |
1.2 Plantejament
Sense perdre de vista els exemples i guies que s’han anat presentant en aquest curs, (penguins o Hawks pel cas de models no supervisats i Titanic pel cas de models supervisats), ara ens trobem amb un conjunt de dades, els terratrèmols o activitat sísmica, amb les quals hem de fer un estudi semblant.
El primer que hem de plantejar és un model no supervisat. Mirant el conjunt de dades, les característiques principals que descriuen els terratrèmols són, la magnitud i la profunditat. No podem fer servir la intensitat ja que això requereix una neteja més exhaustiva que ens portaria a un nou conjunt de dades molt pobre.
El fet que hi hagi tsunami o no, és una variable que depèn més aviat de la situació geogràfica, és a dir, si hi ha aigua a les proximitats del terratrèmol.
Podríem fer servir les variables que descriuen la mortalitat i els danys materials, (tot i que com vam veure a la pràctica 1, depenen d’altres factors, com ara densitat de població, tipus o riquesa del país).
Podem plantejar el problema primer buscant clústers amb un model no
supervisat, per tal de categoritzar els terratrèmols i a partir d’aquí,
fer servir aquesta categorització com a variable objectiu i aplicar un
mètode supervisat.
Alternativament podem buscar primer alguna agrupació i després escollir
una variable objectiu d’entre les que decidim. De fet, seguirem aquesta
darrera opció.
Una primera idea pot ser una classificació dels terratrèmols segons la seva magnitud o la seva profunditat. (Podem discretitzar aquestes variables i mirar si obtenim clústers semblants, tot i que aquesta discretització és més aviat subjectiva i no té perquè coincidir). Una altra possibilitat és escollir com a variable objectiu els danys personals per tal de predir quins terratrèmols causaran aquestes morts.
2 Exercici 1
2.1 Generació d’un model no supervisat.
Les nostres dades són:
str(dades)## 'data.frame': 1137 obs. of 9 variables:
## $ FLAG_TSUNAMI : num 1 1 0 1 0 0 0 0 1 0 ...
## $ YEAR : num 1900 1901 1902 1902 1902 ...
## $ FOCAL_DEPTH : num 25 33 15 33 30 9 100 10 25 20 ...
## $ EQ_PRIMARY : num 8.4 8.2 6.9 7.5 7.7 6.4 7.8 6.2 7.8 6.6 ...
## $ LATITUDE : num 11 40.6 40.7 14 39.9 ...
## $ LONGITUDE : num -66 142.3 48.6 -91 76.2 ...
## $ DEATHS : num 25 18 86 2000 2500 4880 2 4 19000 120 ...
## $ DEATHS_DESCRIPTION: num 1 1 2 4 4 4 1 1 4 3 ...
## $ DAMAGE_DESCRIPTION: num 3 2 4 3 4 4 3 3 4 3 ...És a dir, tenim 1137 observacions i 9 variables.
Primer de tot normalitzem les dades per tal d’evitar biaixos degut a les diferents escale i unitats amb les que traballarem:
# Definim la funció de normalització
nor <- function(x) {(x -min(x))/(max(x)-min(x))}
# Guardem un nou dataset normalitzat
dades <- as.data.frame(lapply(dades, nor))2.2 Comprovació de l’existència de clústers
Abans de buscar el possible nombre de clusters, podem comprovar si el
nostre dataset és apte a tenir clústers.
(Refs. Assessing Clustering Tendency:
- https://www.datanovia.com/en/lessons/assessing-clustering-tendency/
- https://www.r-bloggers.com/2016/10/assessing-clustering-tendency-a-vital-issue-unsupervised-machine-learning/).
Hi ha varies maneres de comprovar-ho. Farem servir el mètode estadístic, l’estadística de Hopkins, que ens dona un valor H que és un indicatiu de l’existència de clústers.
En el nostre cas fem servir la llibreria
factorextra:
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa# Compute Hopkins statistic
res <- get_clust_tendency(dades, n = nrow(dades) - 1, graph = FALSE)
res$hopkins_stat## [1] 0.8763131Un valor de H proper a 1 ens diu que tenim opcions reals d’obtenir clústers, en el nostre cas tenim H = 0.876, per tant podem seguir i esbrinar quants clústers trobem.
2.3 Cerca del nombre de clústers
Desconeixem el nombre de clústers, per tant intentem esbrinar-ho. Comencem amb el mètode del colze, basat en la suma dels quadrats de les distàncies dels punts de cada grup respecte al seu centre (withinss), amb la major separació entre centres de grups (betweenss).
library(cluster)
dissim <- daisy(dades)
resultats_ssw <- rep(0, 10)
resultats_silhouette <- rep(0, 10)
for (i in c(2, 3, 4, 5, 6, 7, 8, 9, 10)) {
set.seed(123)
fit <- kmeans(dades, i)
y_cluster <- fit$cluster
resultats_ssw[i] <- fit$tot.withinss
sk <- silhouette(y_cluster, dissim)
# sk[,1] eś el cluster, sk[,2] és el cluster veí, sk[,3] és l'amplada de la silueta
resultats_silhouette[i] <- mean(sk[, 3])
}
plot(2:10, resultats_ssw[2:10], type = "o", col = "blue", pch = 0, xlab = "Nombre de clústers", ylab = "Suma de quadrats")
Observem que per k = 3 la corba es comença a doblegar,
(tot i que també podem contemplar k = 2 i
k = 4).
Considerem ara el mètode de la silueta que mesura la semblança dels punts d’un clúster respecte a clústers propers. Visualitzem la silueta mitja.
plot(2:10, resultats_silhouette[2:10], type = "o", col = "blue", pch = 0, xlab = "Nombre de clústers", ylab = "Silueta")
Aquesta gràfica no ens diu massa, el que sí veiem que el valor de la silueta és inferior a 0.5.
També es pot fer servir la funció kmeansruns() del
paquet fpc que executarà l’algorisme kmeans com un conjunt de valors i
seleccionarà el valor del número de clústers que millor funcioni d’acord
amb dos criteris: la silueta mitja (asw) i Calinski-Harabasz (“ch”).
library('fpc')
fit_ch <- kmeansruns(dades, krange = 1:10, criterion = "ch")
fit_asw <- kmeansruns(dades, krange = 1:10, criterion = "asw") Podem comprovar el valor amb el que s’ha obtingut el millor resultat
i també mostrar el resultat obtingut per a tots els valors de
k fent servir tots dos criteris.
fit_ch$bestk## [1] 3fit_asw$bestk## [1] 2Si ho visualitzem:
plot(1:10, fit_asw$crit, type = "o", col = "blue", pch = 0, xlab = "Número de clústers",ylab = "Criteri silueta mitja")
Amb el criteri de la silueta mitja obtenim k = 2 però
segons la gràfica k = 3 és proper.
plot(1:10, fit_ch$crit,type = "o", col = "blue", pch = 0, xlab = "Número de clústers", ylab = "Criteri Calinski-Harabasz")
Amb el criteri de Calinski-Harabasz ens dona un valor de 3, però
s’observa que per k = 4 també podria ser un valor
acceptable
Com a resum podem dir que triar el nombre de clústers no és senzill i
menys en el nostre cas que anem una mica a les fosques, no
obtenim un valor concret i s’ha de confrontar amb diferents mètodes. En
el nostre cas no perdrem de vista aquest valor de k = 3 que
hem obtingut amb el mètode les colze i amb el criteri de
Calinski-Harabasz.
2.4 Aplicació de l’algoritme k-means
Tot seguit volem veure com es comporta kmeans. Dels apartats
anteriors hem trobat que el nombre de clústers candidats poden ser
valors de k que van de 2 a 4.
Podem visualitzar inicialment 2, 3, 4 i 5 clústers representant les morts vs la magnitud dels terratrèmols:
set.seed(1)
earthq2clusters <- kmeans(dades, 2)
set.seed(1)
earthq3clusters <- kmeans(dades, 3)
set.seed(1)
earthq4clusters <- kmeans(dades, 4)
set.seed(1)
earthq5clusters <- kmeans(dades, 5)
par(mfrow = c(1,2))
plot(dades[c(4,7)], col = earthq2clusters$cluster, main = "k-means per k = 2", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(4,7)], col = earthq3clusters$cluster, main = "k-means per k = 3", xlab = "Magnitud(EQ_PRIMARY)")
par(mfrow = c(1,2))
plot(dades[c(4,7)], col = earthq4clusters$cluster, main = "k-means per k = 4", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(4,7)], col = earthq5clusters$cluster, main = "k-means per k = 5", xlab = "Magnitud(EQ_PRIMARY)")
Podem observar que la millor agrupació és per k = 2, per
k = 3 també es poden veure agrupacions, però es comencen a
barrejar i per k superior a 3 no obtenim clústers
diferenciats, al contrari els punts es barregen.
Provem altres combinacions amb altres variables.Per exemple morts vs profunditat:
par(mfrow = c(1,2))
plot(dades[c(3,7)], col = earthq2clusters$cluster, main = "Morts vs Profunditat (k = 2)")
plot(dades[c(3,7)], col = earthq3clusters$cluster, main = "Morts vs Profunditat (k = 3)")
Obtenim uns resultats molt semblants als anteriors, per
k = 2 veiem dos clústers, encara que amb alguns punts
sobreposats i per k = 3 els punts estan més barrejats.
Considerem la descripció de morts i danys vs la magnitud i aquest cop
fem la visualització per k = 2, que és la
visualització que farem a partir d’ara, ja que per k = 3
els punt surten barrejats.
par(mfrow = c(1,2))
plot(dades[c(4,8)], col = earthq2clusters$cluster, main = "Descripció morts vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(4,9)], col = earthq2clusters$cluster, main = "Descripció danys vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
En aquest cas es fa difícil apreciar res. Podem justificar-ho dient que tant la descripció de morts com de danys és una divisió poc objectiva, en el sentit que s’ha dividit d’una manera arbitrària i per tant no es presta a una distinció clara.
Continuem amb la visualització de tsunami i profunditat vs magnitud:
par(mfrow = c(1,2))
plot(dades[c(4,1)], col = earthq2clusters$cluster, main = "Tsunami vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(7,1)], col = earthq2clusters$cluster, main = "Tsunami vs Morts")
La darrera gràfica no ens diu res, la primera segueix la dinàmica de
les anteriors. De totes maneres anem veient que la millor separació és
per k = 2.
Si considerem la longitud i la latitud vs la magnitud:
par(mfrow = c(1,2))
plot(dades[c(4,5)], col = earthq2clusters$cluster, main = "Latitud vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(4,6)], col = earthq2clusters$cluster, main = "Longitud vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
Novament no obtenim cap classificació.
Podem inspeccionar la relació amb les morts, que és la variable que potser ens pot funcionar per la classificació.
par(mfrow = c(1,2))
plot(dades[c(7,5)], col = earthq2clusters$cluster, main = "Latitud vs Morts")
plot(dades[c(7,6)], col = earthq2clusters$cluster, main = "Longitud vs Morts")
Podem provar ara amb longitud vs latitud i tsunami vs profunditat:
par(mfrow = c(1,2))
plot(dades[c(4,3)], col = earthq2clusters$cluster, main = "Profunditat vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(3,1)], col = earthq2clusters$cluster, main = "Tsunami vs Profunditat")
2.5 Primeres conclusions
Com a resum, hem vist que el nombre òptim de clústers era
possiblement k = 3, però després de visualitzar varies
combinacions de variables s’observa que per aquest valor de
k els punts surten massa barrejats. En canvi per
k = 2, valor que ens ha donat el criteri de la silueta
mitja, sembla ser el que descriu millor els clústers.
Observem que les variables que més bé mostren els clusters són la magnitud, la profunditat, el nombre de morts i també el fet que hi hagi tsunami o no. Per les altres variables, com longitud i latitud i en especial les descripcions de morts i danys no obtenim resultats apreciables. Una possible explicació dels resultats obtinguts és que les dades de les que disposem han estat sotmeses a una neteja exhaustiva, juntament amb la presència de gran varietat de països on els efectes dels terretrèmols de característiques semblants, poden ser molt diferents. Òbviament una de les altres explicacions, com ja comentat, és la diferència de densitat de població.
3 Exercici 2
3.1 Model no supervisat anterior, però usant una mètrica de distància diferent.
Considerem ara el mateix model que acabem d’estudiar, però usant una mètrica diferent a la mètrica euclidiana que fa servir per defecte l’algoritme kmeans.
La mètrica que farem servir és la distància Manhattan.
library(factoextra)
library(NbClust)
set.seed(123)
res.nb <- NbClust(dades, distance = "manhattan",
min.nc = 2, max.nc = 10,
method = "complete", index ="gap")
res.nb$Best.nc # print the results## Number_clusters Value_Index
## 2.0000 -0.7889fviz_nbclust(dades, kmeans, method = "wss", k.max = 10,
diss = get_dist(dades, method = "manhattan"), nstart = 50)
fviz_nbclust(dades, kmeans, method = "silhouette", k.max = 10,
diss = get_dist(dades, method = "manhattan"), nstart = 50)
Podem observar que el nombre de clusters segueix essent semblant al
que hem determinat, és a dir k = 2 o k = 3
segons les gràfiques anteriors.
Per veure com s’agrupen els clústers fem servir la funció
dist() on indiquem que fem servir la distància Manhattan i
la funció pam(), (Partitioning Around
Medoids).
dist_matrix <- dist(dades, method = "manhattan")
pam_2 <- pam(dist_matrix, k = 2)
pam_3 <- pam(dist_matrix, k = 3)
cluster2 <- pam_2$clustering
cluster3 <- pam_3$clustering
par(mfrow = c(1,2))
plot(dades[c(4,7)], col = cluster2, main = "Morts vs Magnitud (k = 2)", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(4,7)], col = cluster3, main = "Morts vs Magnitud (k = 3)", xlab = "Magnitud(EQ_PRIMARY)")
par(mfrow = c(1,2))
plot(dades[c(3,7)], col = cluster2, main = "Morts vs Profunditat (k = 2)")
plot(dades[c(3,7)], col = cluster3, main = "Morts vs Profunditat (k = 3)")
par(mfrow = c(1,2))
plot(dades[c(4,3)], col = cluster2, main = "Profunditat vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(3,1)], col = cluster2, main = "Tsunami vs Profunditat")
par(mfrow = c(1,2))
plot(dades[c(4,1)], col = cluster2, main = "Tsunami vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
plot(dades[c(7,1)], col = cluster2, main = "Tsunami vs Morts")
par(mfrow = c(1,2))
plot(dades[c(7,5)], col = cluster2, main = "Latitud vs Morts")
plot(dades[c(7,6)], col = cluster2, main = "Longitud vs Morts")
3.2 Comparació dels dos models no supervisats amb mètriques de distància diferents.
Les visualitzacions són molt semblants, hi pot haver alguna lleugera diferència però no és significativa. Fins i tot si ho provem amb la distància del màxim obtenim resultats semblants.
3.3 Conclusions.
Hem fet una nova visualització usant la distància de Manhattan, hem
començat amb valors de k = 2, 3 i hem acabat amb valors de
2 per la mateixa raó que en el cas del primer exercici, els punts es
barrgen com es poden veure a les gràfiques. S’observa que els clústers
són gairebé iguals. Això podria ser degut a la robustesa de les dades o
bé simplement podem dir que les diferències no són significatives al fer
la comparació amb les diferents mètriques. Personalment m’inclinaria per
aquesta darrera, sense entrar en més detall.
4 Exercici 3
4.1 Algorismes DBSCAN i OPTICS
Treballem ara els algoritmes DBSCAN i OPTICS on el paràmetre de
entrada més rellevant és minPts que defineix la mínima
densitat de punts acceptada al voltant d’un centroide i també tenim el
paràmetre eps, (per DBSCAN), que ens indica el radi
respecte els veïns a considerar.
Carreguem la llibreria.
library('dbscan')##
## S'està adjuntant el paquet: 'dbscan'## L'objecte següent està emmascarat per 'package:fpc':
##
## dbscan## L'objecte següent està emmascarat per 'package:stats':
##
## as.dendrogramsummary(dades)## FLAG_TSUNAMI YEAR FOCAL_DEPTH EQ_PRIMARY
## Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.5750 1st Qu.:0.01536 1st Qu.:0.4865
## Median :0.0000 Median :0.7500 Median :0.03379 Median :0.5811
## Mean :0.2067 Mean :0.7048 Mean :0.04859 Mean :0.5785
## 3rd Qu.:0.0000 3rd Qu.:0.8833 3rd Qu.:0.05069 3rd Qu.:0.6757
## Max. :1.0000 Max. :1.0000 Max. :1.00000 Max. :1.0000
## LATITUDE LONGITUDE DEATHS DEATHS_DESCRIPTION
## Min. :0.0000 Min. :0.0000 Min. :0.0000000 Min. :0.0000
## 1st Qu.:0.4930 1st Qu.:0.4981 1st Qu.:0.0000032 1st Qu.:0.0000
## Median :0.7086 Median :0.6463 Median :0.0000222 Median :0.0000
## Mean :0.6409 Mean :0.6075 Mean :0.0059216 Mean :0.1865
## 3rd Qu.:0.7928 3rd Qu.:0.7897 3rd Qu.:0.0001867 3rd Qu.:0.3333
## Max. :1.0000 Max. :1.0000 Max. :1.0000000 Max. :1.0000
## DAMAGE_DESCRIPTION
## Min. :0.0000
## 1st Qu.:0.3333
## Median :0.3333
## Mean :0.5045
## 3rd Qu.:0.6667
## Max. :1.0000Mitjançant les funcions optics i
extractDBSCAN provarem diferents valors per
eps i minPts i ho visualitzarem amb diagrames
d’accessibilitat o reachability plot, (visualització de la
distància d’accessibilitat de cada punt). Les valls representen clústers
(com més profund és la vall, més dens és el clúster), mentre que els
cims indiquen els punts que estan entre les agrupacions.
4.2 Resultats amb
diferents valors d’eps i minPts.
Comencem considerant diferents valors de eps i
minPts = 5:
set.seed(123)
# OPTICS
optics_res = optics(dades, minPts = 5)
# Clústers amb DBSCAN con el valor actual de eps
dbscan_res = extractDBSCAN(optics_res, eps_cl = 0.2)
plot(dbscan_res, main = "minPts = 5, eps = 0.2")
dbscan_res = extractDBSCAN(optics_res, eps_cl = 0.5)
plot(dbscan_res, main = "minPts = 5, eps = 0.5")
dbscan_res = extractDBSCAN(optics_res, eps_cl = 1.2)
plot(dbscan_res, main = "minPts = 5, eps = 1")
Fem el mateix per minPts = 10:
set.seed(123)
# OPTICS
optics_res = optics(dades, minPts = 10)
# Clústers amb DBSCAN con el valor actual de eps
dbscan_res = extractDBSCAN(optics_res, eps_cl = 0.2)
plot(dbscan_res, main = "minPts = 5, eps = 0.2")
dbscan_res = extractDBSCAN(optics_res, eps_cl = 0.5)
plot(dbscan_res, main = "minPts = 5, eps = 0.5")
dbscan_res = extractDBSCAN(optics_res, eps_cl = 1.2)
plot(dbscan_res, main = "minPts = 10, eps = 1.2")
Podem fer-ho també per minPts = 15 i obtenim resultats
semblants.
Donem per bo els valors que ens donen 2 clústers. Els podem
visualitzar representant les gràfiques que hem fet amb kmeans.
Considerem el cas minPts = 10 i eps = 0.5:
set.seed(123)
# OPTICS
optics_res = optics(dades, minPts = 10)
dbscan_res = extractDBSCAN(optics_res, eps_cl = 0.5)
par(mfrow = c(1,2))
hullplot(dades[c(4,7)], dbscan_res, main = "Morts vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
hullplot(dades[c(3,7)], dbscan_res, main = "Morts vs Profunditat")
par(mfrow = c(1,2))
hullplot(dades[c(4,3)], dbscan_res, main = "Profunditat vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
hullplot(dades[c(3,1)], dbscan_res, main = "Tsunami vs Profunditat")
par(mfrow = c(1,2))
hullplot(dades[c(4,1)], dbscan_res, main = "Tsunami vs Magnitud", xlab = "Magnitud(EQ_PRIMARY)")
hullplot(dades[c(7,1)], dbscan_res, main = "Tsunami vs Morts")
par(mfrow = c(1,2))
hullplot(dades[c(7,5)], dbscan_res, main = "Latitud vs Morts")
hullplot(dades[c(7,6)], dbscan_res, main = "Longitud vs Morts")
4.3 Conclusions
Hem vist que per un eps al voltant de 0.5 i fins a
gairebé 1 obtenim 2 clústers, que és el resultat que hem obtingut amb
kmeans. S’observa també que hi ha poques variacions amb
minPts de 5 a 15.
Els resultats de les agrupacions però no són molt correctes ja que hi
ha representacions en les quals s’obtenen clústers superposats. De fet
aquests resultats són molt semblants als que hem trobat amb kmeans, és a
dir, 2 clústers però no ben definits, amb superposicions en els dos
casos estudiats fins ara.
Probablement i com ja hem comentat, per exemple, tenim terratrèmols de
característiques diferents que provoquen els mateixos danys personals i
materials. I de la mateixa manera, tenim terratrèmols amb
característiques iguals que provoquen danys materials i humans molt
diferents.
5 Exercici 4
5.1 Model supervisat i dades
Ara volem plantejar el nostre problema des d’un punt de vista supervisat. El que farem serà escollir una variable objectiu i mirar com es comporta el model amb les dades d’entrenament i test.
La primera cosa que hem de fer és triar la variable objectiu. Podrem pensar amb la magnitud, la profunditat, però ja que disposem de la variable DEATHS_DESCRIPTION, que és una variable discretitzada que conté 4 nivells de danys humans, en decantem per aquesta.
Una altra cosa que podem fer és fer una discretització de la profunditat i la magnitud segons:
- FOCAL_DEPTH (segons la divisió USGS):
- Shallow (superficials): de 0 a 70 km.
- Intermediate (intermedis): de 70 km a 300 km.
- Deep (profunds): de 300 km a 700 km.
#dades_netes$DEPTH <- cut(dades_netes$FOCAL_DEPTH, breaks = c(-1, 70, 300, 700), labels = c("Shallow", "Intermediate", "Deep"))- EQ_PRIMARY (segons la divisió adaptada [UAF] (http://earthquake.alaska.edu/earthquake-magnitude-classes)):
- Minor: fins a magnitud 3.9.
- Light: de 4 a 4.9.
- Moderate: de 5 a 5.9.
- Strong: de 6 a 6.9.
- Great: magnitud 7 o més.
#dades_netes$EQ_MAG_CLASS <- cut(dades_netes$EQ_PRIMARY, breaks = c(0, 4, 5, 6, 7, 10), labels = c("Minor", "Light", "Moderate", "Strong", "Great"))Però després de fer proves i visualitzacions d’arbres ens decantem per deixar-les tal qual, sense discretitzar, el resultats que obtenim fan ús de les característiques pròpies dels terratrèmols.
Primer de tot, eliminem la variable YEAR (considerem que és una
variable que no és significativa de l’activitat sísmica, a més després
de l’exhaustiva neteja realitzada pot donar resultats inesperats).
Tampoc tindrem en compte la variable DEATHS ja que la nostra variable
objectiu és DEATHS_DESCRIPTION.
Un cop dit això, el que farem és crear dos models, un tenint en compte totes les variables i l’altre descartant més variables per tal de simplificar. Amb aquest darrer model seguirem fent els següents exercicis.
Inicialment escollim 7 variables: FLAG_TSUNAMI, FOCAL_DEPTH (profunditat), EQ_PRIMARY (magnitud), LATITUDE, LONGITUDE, DEATHS_DESCRIPTION (variable objectiu) i DAMAGE_DESCRIPTION).
Després descartarem LATITUD i LONGITUD, ens quedarem amb 4 variables per tal d’obtenir un arbre més simplificat.
5.2 Mostres d’entrenament i de prova.
Factoritzem primer les dades que ens interessen:
cols <- c(2, 9, 10, 21, 22, 25, 31)
dades <- as.data.frame(dades_netes[cols])
factor_vars = c("FLAG_TSUNAMI", "DEATHS_DESCRIPTION", "DAMAGE_DESCRIPTION")
dades[factor_vars] <- lapply(dades[factor_vars], factor)Considerem les dades per crear un primer model amb 6 variables i després un altre amb 4 variables, (sense latitud i longitud).
#cols <- c(1, 5, 6, 9, 10, 11)
# Arbre 1
cols <- c(1, 2, 3, 4, 5, 7)
dades_tree_1 <- dades[cols]
# Arbre 2
cols <- c(1, 2, 3, 7)
dades_tree_2 <- dades[cols]Volem dividir el conjunt de dades en un conjunt d’entrenament i un conjunt de prova. El conjunt d’entrenament és el subconjunt del conjunt original de dades utilitzat per a construir un primer model; i el conjunt de prova, el subconjunt del conjunt original de dades utilitzat per a avaluar la qualitat del model.
Dividim aquest subconjunts segons la proporció més usual, 2/3 per al conjunt d’entrenament i 1/3, per al conjunt de prova.
La variable objectiu per la qual classificarem és, com hem dit és DEATHS_DESCRIPTION, (columna 8 del dataset). D’aquesta forma, tindrem un conjunt de dades per a l’entrenament i un per a la validació.
Considerem el primer arbre i definim les dades amb les que treballarem, (afegim seed() per la reproductibilitat).
set.seed(666)
y <- dades[,6]
X_1 <- dades_tree_1De manera dinàmica podem definir una manera de separar les dades en funció d’un paràmetre, en aquest cas del “split_prop”. Definim un paràmetre que controla el split de manera dinàmica en el test.
split_prop <- 3
indexes <- sample(1:nrow(dades), size = floor(((split_prop - 1) / split_prop) * nrow(dades)))
trainX_1 <- X_1[indexes,]
trainy_1 <- y[indexes]
testX_1 <- X_1[-indexes,]
testy_1 <- y[-indexes]Després d’una extracció aleatòria de casos és altament recomanable efectuar una anàlisi de dades mínim per a assegurar-nos de no obtenir classificadors esbiaixats pels valors que conté cada mostra. En aquest cas, verificarem que la proporció de risc és semblant en els dos conjunts:
summary(trainX_1)## FLAG_TSUNAMI FOCAL_DEPTH EQ_PRIMARY LATITUDE
## 0:598 Min. : 0.00 Min. :2.100 Min. :-54.000
## 1:160 1st Qu.: 10.00 1st Qu.:5.700 1st Qu.: 2.471
## Median : 22.00 Median :6.300 Median : 26.474
## Mean : 30.98 Mean :6.365 Mean : 19.493
## 3rd Qu.: 33.00 3rd Qu.:7.100 3rd Qu.: 37.105
## Max. :631.00 Max. :9.500 Max. : 61.017
## LONGITUDE DAMAGE_DESCRIPTION
## Min. :-175.90 1:147
## 1st Qu.: -28.17 2:261
## Median : 51.60 3:169
## Mean : 37.50 4:181
## 3rd Qu.: 103.73
## Max. : 178.29summary(trainy_1)## 1 2 3 4
## 564 44 85 65summary(testX_1)## FLAG_TSUNAMI FOCAL_DEPTH EQ_PRIMARY LATITUDE
## 0:304 Min. : 0.00 Min. :3.500 Min. :-54.000
## 1: 75 1st Qu.: 11.00 1st Qu.:5.800 1st Qu.: 3.011
## Median : 22.00 Median :6.400 Median : 28.700
## Mean : 32.93 Mean :6.413 Mean : 20.158
## 3rd Qu.: 33.00 3rd Qu.:7.100 3rd Qu.: 37.208
## Max. :651.00 Max. :9.100 Max. : 51.613
## LONGITUDE DAMAGE_DESCRIPTION
## Min. :-178.25 1: 63
## 1st Qu.: 13.14 2:140
## Median : 52.70 3: 89
## Mean : 40.09 4: 87
## 3rd Qu.: 102.79
## Max. : 176.51summary(testy_1)## 1 2 3 4
## 269 30 43 37Efectivament s’observa que la proporció és semblant (564/269, 44/30, 85/43, 65/37).
Alternativament, si considerem l’arbre 2, obtenim el següent:
set.seed(666)
X_2 <- dades_tree_2
split_prop <- 3
indexes <- sample(1:nrow(dades), size = floor(((split_prop - 1) / split_prop) * nrow(dades)))
trainX_2 <- X_2[indexes,]
trainy_2 <- y[indexes]
testX_2 <- X_2[-indexes,]
testy_2 <- y[-indexes]
summary(trainX_2)## FLAG_TSUNAMI FOCAL_DEPTH EQ_PRIMARY DAMAGE_DESCRIPTION
## 0:598 Min. : 0.00 Min. :2.100 1:147
## 1:160 1st Qu.: 10.00 1st Qu.:5.700 2:261
## Median : 22.00 Median :6.300 3:169
## Mean : 30.98 Mean :6.365 4:181
## 3rd Qu.: 33.00 3rd Qu.:7.100
## Max. :631.00 Max. :9.500summary(trainy_2)## 1 2 3 4
## 564 44 85 65summary(testX_2)## FLAG_TSUNAMI FOCAL_DEPTH EQ_PRIMARY DAMAGE_DESCRIPTION
## 0:304 Min. : 0.00 Min. :3.500 1: 63
## 1: 75 1st Qu.: 11.00 1st Qu.:5.800 2:140
## Median : 22.00 Median :6.400 3: 89
## Mean : 32.93 Mean :6.413 4: 87
## 3rd Qu.: 33.00 3rd Qu.:7.100
## Max. :651.00 Max. :9.100summary(testy_2)## 1 2 3 4
## 269 30 43 37Obtenim les mateixes proporcions anteriors, (hem triat la mateixa reproductibilitat).
6 Exercici 5
6.1 Creació del model i visualització de l’arbre.
Creem l’arbre de decisió usant les dades d’entrenament (cal no oblidar que la variable outcome és de tipus factor):
library(C50)
library(grid)
library(gridExtra)trainy_1 = as.factor(trainy_1)
model_1 <- C50::C5.0(trainX_1, trainy_1)
plot(model_1, gp = gpar(fontsize = 8))
Com es pot observar, aquest arbre és il·legible. I en canvi ara considerem l’arbre 2 que hem definit més amut i on hem eliminat les variables LATITUD i LONGITUD, tenim el següent:
trainy_2 = as.factor(trainy_2)
model_2 <- C50::C5.0(trainX_2, trainy_2)
plot(model_2, gp = gpar(fontsize = 8))
Aquest darrer arbre s’ha simplificat considerablement. Observem que
s’usen simplement les variables que defineixen intrinsecament els
terratrèmols.
Falta veure ara les característiques i la qualitat.
6.2 Extracció de les regles, comentaris i interpretacions.
Volem analitzar i comentar les dades que podem treure del primer arbre que hem obtingut, un resum del model 1, (en el qual tenim encara les variables de longitud i latitud:
model_1 <- C50::C5.0(trainX_1, trainy_1, rules = TRUE )
summary(model_1)##
## Call:
## C5.0.default(x = trainX_1, y = trainy_1, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jul 30 20:12:00 2024
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 758 cases (7 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (408/34, lift 1.2)
## DAMAGE_DESCRIPTION in {1, 2}
## -> class 1 [0.915]
##
## Rule 2: (711/168, lift 1.0)
## EQ_PRIMARY <= 7.8
## -> class 1 [0.763]
##
## Rule 3: (5/1, lift 12.3)
## FLAG_TSUNAMI = 0
## EQ_PRIMARY > 6.1
## LATITUDE > -7.456
## LATITUDE <= 4.956
## LONGITUDE <= 58.7
## DAMAGE_DESCRIPTION = 3
## -> class 2 [0.714]
##
## Rule 4: (17/11, lift 6.3)
## FLAG_TSUNAMI = 0
## EQ_PRIMARY > 6.1
## LATITUDE <= 28.093
## LONGITUDE <= 58.7
## DAMAGE_DESCRIPTION = 3
## -> class 2 [0.368]
##
## Rule 5: (25/20, lift 3.8)
## EQ_PRIMARY > 8.1
## -> class 2 [0.222]
##
## Rule 6: (41/35, lift 2.8)
## LONGITUDE <= -35.769
## DAMAGE_DESCRIPTION = 4
## -> class 2 [0.163]
##
## Rule 7: (13/2, lift 7.1)
## FLAG_TSUNAMI = 0
## FOCAL_DEPTH > 13
## EQ_PRIMARY > 6.4
## LATITUDE > 28.093
## LONGITUDE > 24
## LONGITUDE <= 122.5
## DAMAGE_DESCRIPTION = 3
## -> class 3 [0.800]
##
## Rule 8: (3, lift 7.1)
## FOCAL_DEPTH > 16
## FOCAL_DEPTH <= 19
## EQ_PRIMARY > 6.8
## LONGITUDE <= 47.448
## DAMAGE_DESCRIPTION = 4
## -> class 3 [0.800]
##
## Rule 9: (3, lift 7.1)
## FOCAL_DEPTH > 16
## FOCAL_DEPTH <= 19
## EQ_PRIMARY > 6.7
## LONGITUDE > 110.2
## DAMAGE_DESCRIPTION = 4
## -> class 3 [0.800]
##
## Rule 10: (11/2, lift 6.9)
## FLAG_TSUNAMI = 0
## FOCAL_DEPTH > 9
## EQ_PRIMARY > 6.1
## LATITUDE > 28.093
## LONGITUDE > 24
## LONGITUDE <= 57.434
## DAMAGE_DESCRIPTION = 3
## -> class 3 [0.769]
##
## Rule 11: (5/1, lift 6.4)
## EQ_PRIMARY > 6.7
## LATITUDE <= -3.1
## LONGITUDE > -35.769
## DAMAGE_DESCRIPTION = 4
## -> class 3 [0.714]
##
## Rule 12: (9/4, lift 4.9)
## EQ_PRIMARY > 7.8
## LATITUDE <= 28.093
## LONGITUDE <= 122.5
## DAMAGE_DESCRIPTION = 3
## -> class 3 [0.545]
##
## Rule 13: (36/23, lift 3.3)
## EQ_PRIMARY > 6.7
## LATITUDE <= 18.339
## DAMAGE_DESCRIPTION = 4
## -> class 3 [0.368]
##
## Rule 14: (9, lift 10.6)
## FOCAL_DEPTH <= 19
## EQ_PRIMARY > 6.8
## LONGITUDE > 47.448
## LONGITUDE <= 110.2
## DAMAGE_DESCRIPTION = 4
## -> class 4 [0.909]
##
## Rule 15: (12/1, lift 10.0)
## FOCAL_DEPTH <= 16
## EQ_PRIMARY > 6.8
## LONGITUDE > -35.769
## DAMAGE_DESCRIPTION = 4
## -> class 4 [0.857]
##
## Rule 16: (26/5, lift 9.2)
## FOCAL_DEPTH > 19
## EQ_PRIMARY > 6.7
## LATITUDE > -3.1
## LONGITUDE > -35.769
## LONGITUDE <= 136.952
## DAMAGE_DESCRIPTION = 4
## -> class 4 [0.786]
##
## Rule 17: (2, lift 8.7)
## FOCAL_DEPTH > 24
## EQ_PRIMARY > 8.1
## LATITUDE <= -24.7
## DAMAGE_DESCRIPTION = 4
## -> class 4 [0.750]
##
## Rule 18: (59/21, lift 7.5)
## EQ_PRIMARY > 6.7
## LONGITUDE > -35.769
## DAMAGE_DESCRIPTION = 4
## -> class 4 [0.639]
##
## Default class: 1
##
##
## Evaluation on training data (758 cases):
##
## Rules
## ----------------
## No Errors
##
## 18 130(17.2%) <<
##
##
## (a) (b) (c) (d) <-classified as
## ---- ---- ---- ----
## 551 4 5 4 (a): class 1
## 36 7 1 (b): class 2
## 51 1 32 1 (c): class 3
## 20 2 5 38 (d): class 4
##
##
## Attribute usage:
##
## 98.68% EQ_PRIMARY
## 72.69% DAMAGE_DESCRIPTION
## 18.87% LONGITUDE
## 13.32% LATITUDE
## 8.71% FOCAL_DEPTH
## 4.49% FLAG_TSUNAMI
##
##
## Time: 0.0 secsObtenim 18 regles, (que es poden llegir en el resum anterior i no reproduïm aquí), i s’observa que el número i percentatge de casos mal classificats en el subconjunt d’entrenament, on hi posa Errors, veiem que es classifiquen erròniament 130 dels 758 casos donats, una taxa d’error del 17.2%.
Si ara fem el mateix per l’arbre 2 (sense lat. i long.):
model_2 <- C50::C5.0(trainX_2, trainy_2, rules = TRUE )
summary(model_2)##
## Call:
## C5.0.default(x = trainX_2, y = trainy_2, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jul 30 20:12:00 2024
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 758 cases (5 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (491/75, lift 1.1)
## EQ_PRIMARY <= 6.7
## -> class 1 [0.846]
##
## Rule 2: (577/89, lift 1.1)
## DAMAGE_DESCRIPTION in {1, 2, 3}
## -> class 1 [0.845]
##
## Rule 3: (127/25, lift 1.1)
## FOCAL_DEPTH > 43
## -> class 1 [0.798]
##
## Rule 4: (73/33, lift 6.4)
## FOCAL_DEPTH <= 43
## EQ_PRIMARY > 6.7
## DAMAGE_DESCRIPTION = 4
## -> class 4 [0.547]
##
## Default class: 1
##
##
## Evaluation on training data (758 cases):
##
## Rules
## ----------------
## No Errors
##
## 4 167(22.0%) <<
##
##
## (a) (b) (c) (d) <-classified as
## ---- ---- ---- ----
## 551 13 (a): class 1
## 40 4 (b): class 2
## 69 16 (c): class 3
## 25 40 (d): class 4
##
##
## Attribute usage:
##
## 85.75% DAMAGE_DESCRIPTION
## 74.41% EQ_PRIMARY
## 26.39% FOCAL_DEPTH
##
##
## Time: 0.0 secsS’observa que el número i percentatge de casos mal classificats en el subconjunt d’entrenament, on hi posa Errors, veiem que es classifiquen erròniament 167 dels 758 casos donats, una taxa d’error del 22.0%.
Les regles que observem són les següents:
- Regla 1:
- Cobreix 491 casos i d’aquests 75 no pertanyen a la classe predita.
- Els terratrèmols que tenen una magnitud menor o igual a 6.7 es classifiquen com a Class = 1. Amb una validesa del 84.6%.
- Regla 2:
- Cobreix 577 casos i d’aquests 89 no pertanyen a la classe predita.
- Els terratrèmols que provoquen uns danys que estàn en els grups 1, 2, 3 es classifiquen com a Class = 1. Amb una validesa del 84.5%.
- Regla 3:
- Cobreix 127 casos i d’aquests 27 no pertanyen a la classe predita.
- Els terratrèmols que tenen profunditat superior a 47 km, es classifiquen com a Class = 1, good. Amb una validesa del 79.8%.
- Regla 4:
- Cobreix 73 casos i d’aquests 33 no pertanyen a la classe predita.
- Els terratrèmols que tenen una profunditat <= 43, una magnitud > 6.7 i uns danys materials = 4, es classifiquen com a Class = 4. Amb una validesa del 54.7%.
6.3 Qualitat del model
Pel primer arbre, com hem dit el número i
percentatge de casos mal classificats en el subconjunt d’entrenament,
veiem que es classifiquen erròniament 130 dels 758 casos donats, una
taxa d’error del 17.2%.
Observem tenim 551 registres ben classificats de la classe 1 i 36, 51 i
20 que han de ser de class 1 i s’han classificat a altres classes. També
tenim 7 registres ben classificats de classe2, 32 de classe 3 i també
tenim 38 registres ben classificats de class 4 i 4, 1 i 1 mal
classificats.
Pel segon arbre, el número i percentatge de casos
mal classificats en el subconjunt d’entrenament, veiem que es
classifiquen erròniament 167 dels 758 casos donats, una taxa d’error del
22.0%.
Observem tenim 551 registres ben classificats de la classe 1 i 40
registres ben classificats de class 4. En canvi no obtenim cap
classificació de les classes 2 i 3.
6.4 Generació la matriu de confusió per a mesurar la capacitat predictiva de l’algoritme, tenint en compte les diferents mètriques associades a aquesta matriu (precisió, sensibilitat, especificitat).
Tot seguit realitzarem una anàlisi de la bondat d’ajust sobre el
conjunt de test amb matriu de confusió, on veurem d’una forma gràfica
com de bé o malament ha funcionat el model.
Un bon model serà el que té valors grans a la diagonal principal i
propers a zero a la resta de posicions de la matriu. Dins la matriu de
confusió, considerem els termes de Sensibilitat (quantifica la
proporció de mostres positives que són classificades com a positives) i
l’Especificitat (quantifica la proporció de mostres negatives
que són classificades com a negatives).
Considerem el primer arbre:
predicted_model_1 <- predict( model_1, testX_1, type="class" )
print(sprintf("La precisió de l'arbre 1 es: %.4f %%", 100 * sum(predicted_model_1 == testy_1) / length(predicted_model_1)))## [1] "La precisió de l'arbre 1 es: 73.3509 %"Per analitzar la matriu de confusió, utilitzarem la funció
CrossTable del paquet gmodels.
library(gmodels)
CrossTable(testy_1, predicted_model_1, prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE, dnn = c('Reality', 'Prediction'))##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 379
##
##
## | Prediction
## Reality | 1 | 2 | 3 | 4 | Row Total |
## -------------|-----------|-----------|-----------|-----------|-----------|
## 1 | 256 | 2 | 5 | 6 | 269 |
## | 0.675 | 0.005 | 0.013 | 0.016 | |
## -------------|-----------|-----------|-----------|-----------|-----------|
## 2 | 23 | 0 | 5 | 2 | 30 |
## | 0.061 | 0.000 | 0.013 | 0.005 | |
## -------------|-----------|-----------|-----------|-----------|-----------|
## 3 | 29 | 1 | 3 | 10 | 43 |
## | 0.077 | 0.003 | 0.008 | 0.026 | |
## -------------|-----------|-----------|-----------|-----------|-----------|
## 4 | 14 | 0 | 4 | 19 | 37 |
## | 0.037 | 0.000 | 0.011 | 0.050 | |
## -------------|-----------|-----------|-----------|-----------|-----------|
## Column Total | 322 | 3 | 17 | 37 | 379 |
## -------------|-----------|-----------|-----------|-----------|-----------|
##
## Per obtenir l’especifitat i la sensibilitat, farem ús de la funció
confusionMatrix() de la llibreria caret.
library(caret)
confusionMatrix(data = predicted_model_1, reference = testy_1, positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2 3 4
## 1 256 23 29 14
## 2 2 0 1 0
## 3 5 5 3 4
## 4 6 2 10 19
##
## Overall Statistics
##
## Accuracy : 0.7335
## 95% CI : (0.686, 0.7774)
## No Information Rate : 0.7098
## P-Value [Acc > NIR] : 0.1682
##
## Kappa : 0.3019
##
## Mcnemar's Test P-Value : 4.639e-08
##
## Statistics by Class:
##
## Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity 0.9517 0.000000 0.069767 0.51351
## Specificity 0.4000 0.991404 0.958333 0.94737
## Pos Pred Value 0.7950 0.000000 0.176471 0.51351
## Neg Pred Value 0.7719 0.920213 0.889503 0.94737
## Prevalence 0.7098 0.079156 0.113456 0.09763
## Detection Rate 0.6755 0.000000 0.007916 0.05013
## Detection Prevalence 0.8496 0.007916 0.044855 0.09763
## Balanced Accuracy 0.6758 0.495702 0.514050 0.73044Fem el mateix pel segon arbre, (sense lat. i long.):
predicted_model_2 <- predict( model_2, testX_2, type="class" )
print(sprintf("La precisió de l'arbre 2 es: %.4f %%", 100 * sum(predicted_model_2 == testy_2) / length(predicted_model_2)))## [1] "La precisió de l'arbre 2 es: 74.1425 %"CrossTable(testy_2, predicted_model_2, prop.chisq = FALSE, prop.c = FALSE, prop.r = FALSE, dnn = c('Reality', 'Prediction'))##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 379
##
##
## | Prediction
## Reality | 1 | 4 | Row Total |
## -------------|-----------|-----------|-----------|
## 1 | 261 | 8 | 269 |
## | 0.689 | 0.021 | |
## -------------|-----------|-----------|-----------|
## 2 | 25 | 5 | 30 |
## | 0.066 | 0.013 | |
## -------------|-----------|-----------|-----------|
## 3 | 32 | 11 | 43 |
## | 0.084 | 0.029 | |
## -------------|-----------|-----------|-----------|
## 4 | 17 | 20 | 37 |
## | 0.045 | 0.053 | |
## -------------|-----------|-----------|-----------|
## Column Total | 335 | 44 | 379 |
## -------------|-----------|-----------|-----------|
##
## L’especifitat i la sensibilitat:
confusionMatrix(data = predicted_model_2, reference = testy_2, positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2 3 4
## 1 261 25 32 17
## 2 0 0 0 0
## 3 0 0 0 0
## 4 8 5 11 20
##
## Overall Statistics
##
## Accuracy : 0.7414
## 95% CI : (0.6942, 0.7848)
## No Information Rate : 0.7098
## P-Value [Acc > NIR] : 0.09553
##
## Kappa : 0.2843
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity 0.9703 0.00000 0.0000 0.54054
## Specificity 0.3273 1.00000 1.0000 0.92982
## Pos Pred Value 0.7791 NaN NaN 0.45455
## Neg Pred Value 0.8182 0.92084 0.8865 0.94925
## Prevalence 0.7098 0.07916 0.1135 0.09763
## Detection Rate 0.6887 0.00000 0.0000 0.05277
## Detection Prevalence 0.8839 0.00000 0.0000 0.11609
## Balanced Accuracy 0.6488 0.50000 0.5000 0.735186.4.1 Comparació i interpretació dels resultats, avaluació la taxa d’error de l’arbre generat, l’eficiència a la classificació
A continuació, analitzem els resultats obtinguts pera cadascun dels models, (dels dos arbres anteriors):
Model 1, (amb les variables de latitud i longitud)
- Train - Errors: classifica erròneament 130 dels 758 casos donats, una taxa d’error del 17.2%.
- Sensibilitat: 95.17% per la Classe 1, 0.0% per la Classe 2, 6.98% per la Classe 3 i 51.35% per la Classe 4.
- Especificitat: 40.00% per la Classe 1, 99.14% per la Classe 2, 95.83% per la Classe 3 i 94.74 per la Classe 4.
- Test - Accuracy: 73.35%
Model 1, (sense les variables de latitud i longitud)
- Train - Errors: classifica erròneament 167 dels 758 casos donats, una taxa d’error del 22.0%.
- Sensibilitat: 97.03% per la Classe 1, 0.0% per la Classe 2, 0.0% per la Classe 3 i 54.05% per la Classe 4.
- Especificitat: 32.73% per la Classe 1, 100% per la Classe 2, 100% per la Classe 3 i 92.98 per la Classe 4.
- Test - Accuracy: 74.14%
6.5 Conclusions.
Com es pot veure, els resultats dels dos models considerats, amb o
sense les variables de LATITUD i LONGITUD, són comparables. Per un
costat la taxa d’error per les dades d’entrenament és lleugerament
inferior quan tenim en comte les variables esmentades.
La precisió per les dades de test és semblant pels dos arbres.
La classificació per la classe 1 és alta, és a dir pels terratrèmols
lleus i sense (o amb molt pocs), morts.
Per les classes 2 i 3 els dos models són pèssims i per la classe 4 està
al voltant del 50%, que ho qualifiquem de dolent.
7 Exercici 6
7.0.1 Es prova amb una variació o un altre enfocament algorítmic.
Boosting:
Considerem ara l’adaptative boosting que és una altra característica potent incorporada a C5.0, està basat en el treball de Rob Schapire i Yoav Freund. La idea és generar diversos classificadors (ja siguin arbres de decisió o conjunts de regles) en lloc d’un sol. Cada cop que s’ha de classificar un cas nou, cada classificador vota per la seva classe prevista i els vots es compten per determinar la classe final. (https://www.rulequest.com/see5-unix.html#RULES).
Nota: ho farem només per el primer dels models anteriors, el qual hi tenim les variables de latitud i longitud.
model2 <- C50::C5.0(trainX_1, trainy_1, trials = 45)
#plot(model2, gp = gpar(fontsize = 8))Novament obtenim un arbre que no és llegible per la seva densitat, per tant no el mostrem. El que sí mostrem són les noves prediccions d’aquest nou arbre.
predicted_model2 <- predict( model2, testX_1, type="class" )
print(sprintf("La precisió de l'arbre és: %.4f %%",
100 * sum(predicted_model2 == testy_1) / length(predicted_model2)))## [1] "La precisió de l'arbre és: 72.8232 %"Podem veure que la precisió de l’arbre gairebé és semblant, (ha passat de 73.35% a 72.82%).
Tot seguit inspeccionem la matriu de confusió.
# Obtenció de la sensibilitat i especificitat
confusionMatrix(data = predicted_model2, reference = testy_1,
positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2 3 4
## 1 252 21 28 11
## 2 4 0 0 2
## 3 9 5 5 5
## 4 4 4 10 19
##
## Overall Statistics
##
## Accuracy : 0.7282
## 95% CI : (0.6805, 0.7724)
## No Information Rate : 0.7098
## P-Value [Acc > NIR] : 0.2321
##
## Kappa : 0.3167
##
## Mcnemar's Test P-Value : 1.693e-05
##
## Statistics by Class:
##
## Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity 0.9368 0.00000 0.11628 0.51351
## Specificity 0.4545 0.98281 0.94345 0.94737
## Pos Pred Value 0.8077 0.00000 0.20833 0.51351
## Neg Pred Value 0.7463 0.91957 0.89296 0.94737
## Prevalence 0.7098 0.07916 0.11346 0.09763
## Detection Rate 0.6649 0.00000 0.01319 0.05013
## Detection Prevalence 0.8232 0.01583 0.06332 0.09763
## Balanced Accuracy 0.6957 0.49140 0.52987 0.73044S’observa que ha disminuït la sensibilitat respecte la classe 1, però ha augmentat llugerament la de la classe 3, (11.63%).
Cart
Considerem ara arbres de decisions generats per l’algoritme CART
(Classification And Regression Trees) mitjançant la funció
rpart().
Aquesta funció implementa arbres de classificació utilitzant un punt de vista de partició recursiva que divideix el conjunt de dades en subconjunts homogenis, basats en els valors de les variables predictores. Aquest subconjunts es representen amb nodes en un arbre, on els nodes terminals contenen les prediccions basades en la majoria o el promig de les observacions d’aquest subconjunt.
Les regles de decisió generades per el model CART, generalment es visualitzen com un arbre binari.
if(!require(rpart)){
install.packages('rpart', repos='http://cran.us.r-project.org')
library(rpart)
}## S'està carregant el paquet requerit: rpartif (!require('rpart.plot')) install.packages('rpart.plot')## S'està carregant el paquet requerit: rpart.plotlibrary('rpart.plot')
set.seed(123)
# Model rpart
model_rpart <- rpart(dades$DEATHS_DESCRIPTION ~ FLAG_TSUNAMI + FOCAL_DEPTH + EQ_PRIMARY +
LATITUDE + LONGITUDE + DAMAGE_DESCRIPTION,
data = dades_tree_1,
method = "class")
# Plot rpart
rpart.plot(model_rpart, type = 4, extra = 102)
set.seed(123)
# Model rpart
predicted_rpart <- predict(model_rpart, testX_1, type = "class")
print(sprintf("Model rpart: La precisió de l'arbre és: %.4f %%",
100 * sum(predicted_rpart == testy_2) / length(predicted_rpart)))## [1] "Model rpart: La precisió de l'arbre és: 74.6702 %"# Obtenció de la sensibilitat i especificitat
confusionMatrix(data = predicted_rpart, reference = testy_1,
positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 2 3 4
## 1 262 28 33 16
## 2 0 0 0 0
## 3 0 0 0 0
## 4 7 2 10 21
##
## Overall Statistics
##
## Accuracy : 0.7467
## 95% CI : (0.6998, 0.7897)
## No Information Rate : 0.7098
## P-Value [Acc > NIR] : 0.06186
##
## Kappa : 0.2862
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity 0.9740 0.00000 0.0000 0.56757
## Specificity 0.3000 1.00000 1.0000 0.94444
## Pos Pred Value 0.7729 NaN NaN 0.52500
## Neg Pred Value 0.8250 0.92084 0.8865 0.95280
## Prevalence 0.7098 0.07916 0.1135 0.09763
## Detection Rate 0.6913 0.00000 0.0000 0.05541
## Detection Prevalence 0.8945 0.00000 0.0000 0.10554
## Balanced Accuracy 0.6370 0.50000 0.5000 0.75601En aquest cas, la precisió de l’arbre ha augmental lleugerament, de 73.35% a 74.67%.
En quant a la sensibilitat, ha augmentat per la classe 1 i la classe 4. Però per les classes 2 i 3 seguim amb classificació zero.
8 Exercici 7
8.1 Possibles limitacions tenen les dades que has seleccionat per obtenir conclusions amb els models (supervisat i no supervisat)
Hem fet un estudi dels terratrèmols arreu del món considerant un conjunt de dades altament incomplet, degut bàsicament a que els registres que conté el dataset van del 2150 BC fins a l’actualitat i les dades anteriors a l’existència d’aperells d’enregistrament d’activitat sísmica s’han descartat. De fet s’han descartat les dades anteriors al 1900 i algunes de les posteriors per incompletesa. Aquesta neteja exhaustiva ha pogut influir als resultats que hem obtingut.
Hi ha una característica comuna que hem obtingut en els dos models estudiats, el no supervisat i el supervisat.
A la primera part, model no supervisat hem obtingut els millors
resultats amb 2 clústers. En la segona part hem considerat una variable
objectiu (DEATHS_DESCRIPTION) que engloba 4 graus de danys humans. Hem
obtingut una classificació dels terratrèmols segons aquesta variable, en
la qual bàsicament només diferenciava la Classe 1 amb bons resultats i
la Classe 4 amb més mal resultat.
En canvi per la classe 2 i 3 no ens ha funcionat. És a dir, el nostre
model no supervisat ha identificat 2 clústers i el model supervisat ha
identificat 2 classes.
8.2 Possibles riscos de fer servir el model.
És evident que el model no funciona, pels resultats obtinguts. Ara bé, el fet de la falta de classificació de les classes 2 i 3 és un indicatiu de que primer, la divisió de danys humans potser no és la més correcta en el sentit que potser es podria fer una classificació més ample entre les classe 2, 3 i 4.
Per altra part tenim la falta d’informació de la densitat de població. Aquest fet pot introduir un soroll que afecta a la classificació ja que per exemple, terratrèmols que potser es classificarien com a classe 4 van a parar a altres classes inferiors.
Un altre fet, des del meu punt de vista, clau en la mala classificació dels grups 2 i 3, és que tenim una barreja de països de diferents nivells de preparació o prevenció de terratrèmols i per tant amb conseqüències molt variades. És a dir, tenim terratrèmols amb les mateixes característiques intrínseques però amb clasificació de danys humans molt variat. De la mateixa manera, tenim terratrèmols amb característiques molt diferents, (e.g. magnitud), amb danys humans molt semblant, que correspondrien a la mateixa classe.
Això ens porta finalment a preguntar-nos si aquest model podria funcionar si, en comptes de considerar terratrèmols globals, consideressim terratrèmols per diferents països o regions. En aquest cas potser es minimitzarien les diferèncias o les influències que fan que tinguem aquesta mala classificació en el nostre cas estudiat.